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S+ WAVELETS:  Algorithms  and  Technical 

Details 

Abstract 


A  complete  description  is  given  for  the  algorithms  in  S+Wavelets  software 
toolkit  for  wavelet  and  cosine  packet  analysis.  These  algorithms  include  wavelet 
transforms,  wavelet  packet  transforms,  cosine  packet  transforms,  and  non-decimated 
wavelet  transforms.  Implementations  for  the  transforms  and  their  inverses  are  given 
for  a  variety  of  boundary  treatment  rules,  including  periodic,  reflection,  interval 
wavelets  (Cohen  et  al.  [CDV93]),  and  zero/polynomial  extension.  In  addition, 
modifications  to  the  standard  algorithms  are  given  to  handle  signals  or  images  with 
dimensions  not  divisible  by  a  power  of  two. 


Keywords:  Discrete  wavelet  transform,  discrete  cosine  transform,  wavelet  packet 
transform,  boundary  rules,  algorithms,  S+WAVELETS. 
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1  Introduction 

This  paper  describes  the  algorithms  in  the  S+Wavelets  toolkit  for  the  wavelet  and 
cosine  packet  analysis.  S+Wavelets  is  an  S-PLUS  module  for  wavelet  analysis  of 
signals  and  images.  See  Bruce  and  Gao  [BG94,  BG95]  for  details  concerning  the 
S+Wavelets. 

The  classical  one-dimensional  discrete  wavelet  transform  maps  a  continuous  sig¬ 
nal  f{t)  €  L2{R)  to  the  wavelet  coefficients  {sjr,n  :  n  €  Z}  and  {dj^n  '•  j  <  J,n  £  Z}. 
The  signal  is  defined  on  the  entire  real  line  R  and  the  transform  has  an  infinite 
number  of  coefficients.  Given  coefficients  {sjk,n  :  G  Z},  it  is  straightforward  to 
compute  the  coefficients  at  successively  coarser  levels  j  >  k  using  the  pyramid  algo¬ 
rithm  [Mal89]. 

Application  of  the  discrete  wavelet  transform  to  problems  involving  real  data 
involves  several  conceptual  and  computational  hurdles.  Typically,  we  are  only  con¬ 
cerned  with  finitely  sampled  signals  /i,  /2,  •  •  •,  /n-  The  length  of  the  signal  sequence 
N  is  not  necessarily  a  power  of  two,  and  we  may  need  to  keep  all  of  the  data  for  the 
analysis.  In  this  paper,  we  address  these  issues  by  discussing  a  suite  of  algorithms 
for  applying  the  discrete  wavelet  transform  and  the  wavelet  packet  transform  to 
finite  length  signals  of  arbitrary  length. 

To  handle  finite  signals,  special  treatment  is  required  at  the  boundaries.  In 
S-h Wavelets,  several  boundary  treatment  rules  are  available.  These  rules  fall  into 
one  of  three  categories: 

1.  Infinite  Extension  Model:  The  finite  signal  is  extended  to  an  infinite  signal 
•  •  •,  /_i,  /o,  fi,  ■  ■  ■  and  the  the  classical  wavelet  transform  is  applied  to  the 
infinite  sequence  of  coefficients.  The  transform  coefficients  are  then  selected 
from  the  infinite  transform  which  correspond  to  the  original  finite  signal.  Infi¬ 
nite  extension  boundary  rules  in  S-|- WAVELETS  are  periodic,  reflection,  zero, 
and  polynomial  extension. 

2.  Recursive  Extension:  An  ad  hoc  boundary  rule  can  be  defined  by  recursively 
extending  the  signal  at  each  filtering  step  in  the  wavelet  or  wavelet  packet 
transform.  Infinite  extension  requires  either  storage  of  extra  coefficients  or 
extremely  messy  bookkeeping.  As  a  simple  alternative  to  infinite  extensive, 

,,  S-|- Wavelets  offers  zero  and  polynomial  recursive  extension  rules. 

3.  Wavelets  on  an  Interval  Model:  Cohen,  Daubechies,  and  Vial  [CDV93] 
formally  defined  a  new  wavelet  transform  on  the  compact  set  [a,  6].  This 
transform  may  be  different  from  the  classical  wavelet  transform  defined  on  the 
real  line  TZ.  The  periodic  and  reflection  infinite  extension  rules  can  also  be 
viewed  as  wavelet  transforms  restricted  to  a  compact  interval. 
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Cosine  packet  analysis  is  commonly  known  as  local  cosine  analysis,  and  was 
first  introduced  by  Ronald  Coifman  and  Yves  Meyer  [CM91].  The  term  “cosine 
packets”  was  coined  by  David  Donoho,  another  leading  wavelets  researcher,  because 
cosine  packet  analysis  is  a  mirror  image  of  wavelet  packet  analysis  [CMQW90]. 
The  difference  is  that  localized  cosine  functions  are  used  instead  of  wavelet  packet 
functions. 

The  Fourier  cosine  transform  (FCT)  of  a  signal  f{t)  is  given  by 


(1) 


The  discrete  cosine  transform  is  a  discretized  version  of  (1).  There  are  four  com¬ 
monly  used  orthogonal  discrete  cosine  transforms:  DCT-I,  DCT-II,  DCT-III,  and 
DCT-IV.  S-b Wavelets  provides  functions  for  two  of  these  transforms:  DCT-II 
and  DCT-IV. 

Cosine  packet  analysis  is  a  generalization  of  the  discrete  cosine  transform  (DCT). 
The  DCT  is  widely  used  in  signal  processing  and  image  processing,  and  is  particu¬ 
larly  valuable  for  coding  and  data  compression  applications.  One  drawback  of  the 
DCT  is  the  abrupt  “cutoff”  implicit  in  dividing  the  signal  into  disjoint  blocks.  This 
can  cause  undesirable  block  effects,  such  as  “Gibbs  phenomena.”  To  avoid  the  prob¬ 
lems  caused  by  the  abrupt  cutoff,  Coifman  and  Meyer  [CM91]  introduced  a  new  type 
of  localized  cosine  transform  with  smooth  cutoffs  (tapers). 

A  cosine  packet  function  is  obtained  by  damping  a  cosine  function  down  to  zero 
on  an  interval  I  using  a  taper  function  or  bell  function  Type  II  and  type  IV 
cosine  packet  functions  with  frequency  k  defined  on  the  interval  I  =  [a,  6]  are  given 

by 
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where  A/  =  6  —  a.  In  order  to  define  orthogonal  transforms,  a  very  special  kind  of 
tapering  function  is  needed. 

'Section  2  gives  the  “vanilla”  algorithm  for  computing  the  discrete  wavelet  trans¬ 
form  and  wavelet  packet  transform.  Treatment  of  the  boundaries  is  ignored  in  the 
vanilla  algorithm.  Section  3  gives  an  overview  of  the  boundary  rules  for  the  dis¬ 
crete  wavelet  transform  and  for  the  wavelet  packet  transform.  Section  4  gives  the 
wavelet  convolution  and  down-sampling  algorithms.  Section  5  gives  the  wavelet 
convolution  and  up-sampling  algorithms.  Section  6  gives  the  algorithm  for  comput¬ 
ing  non-decimated  (over-sampled)  wavelet  transform,  along  with  an  algorithm  for 
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wavelet  convolution  and  dilation.  Section  7  describes  the  algorithms  used  to  com¬ 
pute  the  discrete  cosine  transforms  (DCT-II  and  DCT-IV).  The  Section  7.3  describes 
the  boundary  extension  rules  for  cosine  packets.  Section  8  describes  the  algorithm 
for  computing  a  cosine  packet  transforms.  The  inverse  cosine  packet  transform  is 
discussed  in  section  9.  The  equations  for  the  tapers  are  given  in  Appendix  section  B. 

This  paper  assumes  the  reader  is  familiar  with  the  beisic  methods  and  algo¬ 
rithms  for  wavelets  and  wavelet  packets.  For  background  on  wavelets,  refer  to 
[Str89,  Dau92,  Chu92,  JLS94].  For  background  on  wavelet  packets,  refer  to  [CW92, 
CMW92,  Wic94].  For  further  discussion  regarding  the  algorithms  ajid  methods  be¬ 
hind  cosine  packet  analysis  (local  cosine  analysis),  refer  to  the  book  by  Wickerhauser 
[Wic94]. 

2  Wavelet  and  Wavelet  Packet  Transforms:  Vanilla 
Algorithms 

The  discrete  wavelet  transform  (DWT)  and  the  inverse  discrete  wavelet  transform 
(IDWT)  use  Mallat’s[Mal89]  remarkable  fast  pyramid  algorithms.  The  wavelet 
packet  transform  (WPT)  and  inverse  wavelet  packet  transform  (IWPT)  involve  a 
generalization  of  the  pyramid  algorithm  [CMQW90].  Mallat’s  wavelet  pyramid  algo¬ 
rithm  has  its  roots  in  the  pyramid  algorithm  of  Burt  and  Adelson  [BA83].  Chapter 

3  of  [Mey93]  gives  a  historical  perspective  of  the  pyramid  algorithm. 

The  forward  algorithms  involve  convolution  with  low-pass  and  high-pass  analysis 
filters  followed  by  down-sampling  (decimation)  operator.  The  backwards  (inverse) 
algorithms  involve  convolution  with  low-pass  and  high-pass  synthesis  filters  followed 
by  up-sampling  (zero-padding)  operator.  The  convolution  and  down-sampling  oper¬ 
ators  are  described  in  section  4,  and  the  convolution  and  up-sampling  operators  are 
described  in  section  5.  In  this  section,  we  describe  the  “vanilla”  wavelet  transform 
ajid  wavelet  packet  transform  algorithms. 


2.1  Forward  Algorithm  to  Compute  the  DWT 

The  DWT  pyramid  algorithm  is  represented  by  figure  1.  The  basic  components  in 
each  stage  of  the  pyramid  are  two  analysis  filters — a  low-pass  filter  L  and  a  high- 
pass  filter  H — and  a  decimation-by-two  operation.  The  down-sampling  (decimation) 
operation,  indicated  by  a  downward  pointing  arrow  with  the  number  2,  consists  of 
deleting  every  other  value  of  the  filter  outputs. 

Start  with  a  signal  z  =  {zi,Z2,. . Zn)'  ■  Let  Wp^i  be  the  convolution  and  down- 
sampling  operator  for  filter  F.  The  pyramid  algorithm  of  J  multiresolution  levels 
consists  of  the  following  steps; 
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Figure  1:  DWT  pyramid  algorithm. 


[0]  Initialize  the  input  So  =  (5o,i,>So,2j •  •  •  ,-so,n)'  to  the  pyramid  algorithm  as  fol¬ 
lows: 

SQ,i  =  Zi  i  =  l,2,---,n. 

Initialize  the  level  index  i  =  1. 

[1]  Apply  the  convolution  and  down-sampling  operator  to  sy-i  with  the  high  pass 
filter  H  to  obtain  the  level  j  detail  coefficients 


d,  = 


Store  the  dj  coefficients. 

[2]  Apply  the  convolution  and  down-sampling  operator  to  Sj-i  with  the  low  pass 
filter  L  to  obtain  the  level  j  smooth  coefficients 

Sj  = 

If  j  <  J,  then  increment  j  =  j  +  1^  and  go  to  step  1.  Otherwise,  store  the  Sj 
coefficients  and  quit. 

The  output  of  the  algorithm  is  the  set  of  DWT  coefficients  di,  d2,  •  •  •,  dj,  and  sj. 

2.2  Backwards  Algorithm  to  Compute  the  IDWT 

The  inverse  discrete  wavelet  transform  (IDWT)  algorithm,  represented  by  figure  2, 
inverts  the  DWT  pyramid  algorithm  in  a  straightforwaxd  manner.  The  basic  compo¬ 
nents  in  each  stage  of  the  reconstruction  are  the  synthesis  low  and  high  pass  filters 
L’  and  H“  and  an  up-sample-by-two  operation.  The  up-sample  operation,  indicated 
by  an  upward  pointing  arrow  with  the  number  2,  consists  of  inserting  zeros  between 
every  other  value  of  the  filter  inputs. 

Start  with  the  DWT  coefficients  di,  d2,  •  ■  dj,  and  sj.  Let  be  the  con¬ 
volution  and  up-sampling  operator  for  filter  F.  The  backwards  pyramid  algorithm 
proceeds  as  follows: 

[0]  Initialize  the  level  index  j  =  J. 
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Figure  2:  Reconstruction  algorithm  for  the  discrete  wavelet  transform. 

[1]  Apply  the  convolution  ajid  up-sampling  operators  to  the  level  j  coefficients  to 
obtain  the  level  j  —  1  smooth  coefficients: 

Si-,  =  Wuisi)  +  w„,,{d,). 

[2]  If  ;■  >  1,  then  go  to  step  1.  Otherwise,  for  i  =  1,  set  the  output  signal  z  to 

Z,’  =  So,i  Z  =  1,  2,  •  •  •  ,  TZ. 


2.3  Wavelet  Packet  Table 

The  pyramid  algorithm  for  wavelet  packet  table  is  represented  by  figure  3.  The 
basic  components  in  each  stage  of  the  pyramid,  similar  to  DWT  case,  are  two  anal¬ 
ysis  filters — a  low-pass  filter  L  and  a  high-pass  filter  H — and  a  decimation-by-two 
operation.  The  down-sampling  (decimation)  operation,  indicated  by  a  downward 


wO.O 


w1.0 


w2.0 


w2.1 


w1.1 


W2.2 


w2.3 


'  Figure  3:  Pyramid  algorithm  for  wavelet  packet  table. 

pointing  arrow  with  the  number  2,  consists  of  deleting  every  other  value  of  the  filter 
outputs. 

Start  with  a  signal  z  =  (zi,  Z2, . . . ,  Zn)' .  Let  Wf,[  be  the  convolution  and  down- 
sampling  operator  for  filter  F.  The  pyramid  algorithm  of  J  multiresolution  levels 
consists  of  the  following  steps: 
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[0]  Initialize  the  input  Wo,o  =  {wo,o,i,wo,o,2^  •  •  •  i^o.o.n)^  to  the  pyramid  algorithm 
35  follows: 

WQ,o,i  =  2i  i  =  l,2,---,n. 

Initialize  the  level  index  j  =  1. 

[1]  For  each  b  =  0,1,..., 2^“^  —  1,  apply  the  convolution  and  down-sampling 
operator  to  with  the  high  pass  filter  H  and  the  low  pass  filter  L  to 

obtain  the  level  j  coefficients.  The  coefficients  axe  stored  in  sequency  order:  if 
h  is  even. 


and  if  b  is  odd. 


Wj,26  =  FFi.iCWi-l,!,) 

Wj.2i+i  = 


Wj,26  = 

Wj,2i+1  =  WL,x(Wj_i,6). 


[2]  If  j  <  J,  then  increment  j  =  j  +  1,  and  go  to  step  1.  Otherwise  quit. 

The  output  of  the  algorithm  is  the  set  of  wavelet  packet  coefficients  which  are 

stored  as  a  long  vector  (except  when  boundary= '  'infinite' ')  .  ,wj^2-’-i)- 


2.4  Inverse  Wavelet  Packet  Transform 

The  inverse  wavelet  packet  transform  (IWPT)  algorithm  inverts  the  wavelet  packet 
pyramid  algorithm.  Start  with  a  table  of  wavelet  packet  coefficients  {wj,;,}  where 
j  €  {0, . . . ,  J}  and  6  G  {0, . . . ,  2-'  —  1}  (to  handle  matching  pursuits  these  do  not 
necessarily  have  to  correspond  to  an  orthogonal  basis). 

[0]  Initialize  a  logical  vector  of  length  2^  J  -|- 1)  —  1  indicating  whether  the  (j,  6) 
block  is  in  the  wavelet  packet  table.  Initialize  the  current  level  j  =  J . 

[1]  For  each  6  =  0, 1, ... ,  —  1,  if  block  {j,  b)  is  in  the  table,  apply  the  convo- 

,  lution  and  up-sampling  operator  using  the  appropriate  filter.  If  6%%4  =  0  or 

6%%4  =  3,  where  %%  is  the  modulus  operator,  then  used  a  low-pass  filter  L. 
Otherwise  use  a  high-pass  filter  H.  Add  the  result  to  block  [j  —  1,  [V^J)  and 
mark  that  block  as  in  the  wavelet  packet  table. 

[2]  If  j  >  0,  then  decrement  j  =  j  —  1,  and  go  to  step  1.  Otherwise  quit. 

The  output  of  the  algorithm  is  the  reconstructed  vector. 
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3  Overview  of  Boundary  Rules  for  Wavelets 

There  axe  five  basic  types  of  boundary  treatment  rules  for  wavelet  analysis  supported 
in  S+Wavelets:  periodic,  reflection,  recursive  extension  (zero  and  polynomial), 
interval,  and  infinite  extension  (zero  and  polynomial).  These  are  described  below. 

Periodic  Extension 

The  original  series  zx^zi-,. . .  ,2„  is  assumed  to  be  n  periodic,  so  Zi  =  2,®n-  Equiva¬ 
lently,  the  wavelets  are  assumed  to  be  periodic  on  the  interval  [0,  n]. 

Reflection  Extension 

The  original  series  z\,zi,. ..  ,Zn  is  reflected  at  the  boundaries  and  then  periodically 
extended  using  the  algorithm  given  by  [Bri92].  The  reflection  boundary  correc¬ 
tion  gives  perfect  reconstruction  only  for  symmetric  wavelets  (i.e.,  the  biorthogonal 
wavelets). 

Recursive  Extension 

At  each  filtering  step,  the  coefficients  are  padded  at  the  beginning  and  the  end  of  the 
signal  with  zeros  or  with  a  polynomial  extension.  Three  polynomial  extension  rules 
are  implemented  in  S-f- WAVELETS:  polyO  (the  first  and  last  value  are  repeated), 
polyl  (the  beginning  and  the  end  of  the  signal  are  extended  using  a  polynomial 
of  degree  one  fit  to  the  first  two  and  last  two  values  respectively),  and  poly2  (the 
beginning  and  the  end  of  the  signal  are  extended  using  a  polynomial  of  degree  two 
fit  to  the  first  three  and  last  three  values  respectively). 

Interval  Wavelets 

This  rule  corresponds  to  the  special  wavelet  functions  at  the  boundaries  defined  by 
Cohen,  Daubechies  and  Vial[CDV93].  The  boundary  wavelets  are  zero  outside  of 
the  range  of  the  data.  The  transform  retains  the  orthogonality  properties  of  the 
“classical”  wavelet  transform  and  is  numerically  stable. 

For  the  interval  wavelet,  then  a  preconditioning  transform  can  be  (optionally) 
apjilied  to  the  signal  before  applying  the  interval  wavelet  transform.  The  precon¬ 
ditioning  transform  preserves  the  “vanishing  moment”  property  of  wavelets  at  the 
expense  of  introducing  an  additional  non-orthogonal  transform.  See  [CDV93]  for 
details. 
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Infinite  Extension 

The  original  series  zi,Z2i...  ,Zn  is  extended  once  at  the  beginning  of  the  filtering 
procedure  using  a  zero  or  polynomial  extension.  This  is  similar  to  the  zero  and 
polynomial  rules  above,  for  which  the  coefficients  are  extended  at  each  filter  step. 
For  the  polynomial  extension,  a  polynomial  of  degree  pdeg  is  fit  using  a  fraction 
pfrac  of  the  data  (0  <  pfrac  <  1). 

The  infinite  boundary  rule  produces  an  infinite  set  of  wavelet  coefficients.  How¬ 
ever,  only  a  finite  number  of  coefficients  is  stored.  Since  a  polynomial  extension  rule 
is  used,  the  stored  coefficients  can  be  used  to  compute  the  remaining  coefficients. 

Unlike  the  other  boundary  rules,  the  infinite  boundary  treatment  is  an  “expan¬ 
sionist”  transform.  This  means  that  you  end  up  with  more  coefficients  than  original 
sample  values.  For  a  series  of  length  n,  you  obtain  n-\-p  wavelet  coefficients  where 
0  <  p  ^  n,  independent  of  n.  As  a  result,  the  output  data  structure  for  a  transform 
computed  with  the  infinite  boundary  rule  is  different  than  transforms  computed  with 
other  rules  (the  transform  is  a  “crystal  list”  object  rather  than  a  “crystal  vector” 
object:  see  [BG95]). 


4  Convolution  and  Downsampling  Algorithms 

Let  X  =  (xi, . . . ,  x„)  be  the  input  signal  and  let  /  =  (/i, . . . ,  f^)  be  the  filter.  The 
“generic”  convolution  and  down-sampling  operation  is  given  by 

m 

yk  =  Yl  fi^2k+i-2  (4) 

»=1 


Since  Xj  is  only  defined  for  i  =  1, . . . ,  ra,  the  summation  in  (4)  needs  to  be  modified 
at  the  beginning  and  end  of  the  signal.  The  way  in  which  the  summation  is  handled 
at  the  boundaries,  as  well  as  the  length  of  the  output  signal,  depends  on  the  bound¬ 
ary  rule.  This  section  gives  the  algorithms  for  the  convolution  and  down-samphng 
operator  (4)  for  the  boundary  rules  of  section  3. 

Define  the  operator  “convdown.general”  by  (4)  with  k  restricted  to  values  for 
which  the  filter  does  not  extend  beyond  the  ends  of  the  input  signal 


k  =  1,2,.. 


n  —  m  +  2 
2 


The  length  of  the  output  signal  of  the  convdown.general  operator,  £  =  , 

is  in  general  smaller  than  the  desired  length.  In  order  to  obtain  enough  output 
coefficients  ([n/2J  or  [n/2j  +  1),  extra  values  (roughly  m  —  2)  must  be  padded  to 
X  before  calling  the  convdown.general  operator. 

The  extra  values  can  be  all  padded  at  the  begining  or  at  the  end  of  X,  or  partly 
at  the  begining  and  partly  at  the  end.  In  S-h WAVELETS,  if  the  number  of  padding 
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values  is  even,  an  equal  number  of  values  are  padded  at  the  beginning  and  end.  If 
the  number  of  padding  values  is  odd,  then  one  more  value  is  padded  at  the  end. 

In  the  algorithms  below  (except  for  the  reflection  boundary  rule)  if  /  is  of  odd 
length,  we  pad  a  zero  to  /  to  make  it  even  length.  The  filter  padding  rule  is  as 
follows:  if  /  is  a  low-pass  filter,  /  =  (/,  0);  if  /  is  a  high-pass  filter,  /  =  (0,/).  For 
the  reflection  case,  special  filtering  rules  apply  to  odd  length  filters. 

4.1  Periodic 

The  input  signal  X  is  assumed  to  be  periodic  with  period  n.  The  algorithm  consists 
of  the  following  two  steps: 

[1]  extend  the  input  signal 

X  (x,j— p.^1 ,  .  .  .  ,  X,i,  ■  *  *  T  ^1)  •  *  *  j  ^p)' 

[2]  apply  the  general  convdown  operator: 

Y  =  convdown.general(X,  /). 

Note:  If  X  is  assumed  to  be  periodic,  then  Y  is  also  periodic.  When  n  is  even, 
the  period  of  Y  is  n/2.  When  n  is  odd,  the  period  of  Y  is  still  n.  Hence,  sample  size 
of  the  original  signal  is  restricted  to  multiples  of  2"^  for  the  periodic  boundary  rule. 


4.2  Reflection 


The  input  signal  X  is  reflected  and  then  (implicitly)  periodized  using  the  approach 
described  by  [Bri92].  When  a  symmetric/antisymmetric  filter  is  applied,  the  filtered 
vector  Y  is  of  the  same  reflection/periodicity  property.  Following  the  notation  of 
[Bri92],  we  use  the  “£'(1,1)  extension”  for  odd  length  filter  (i.e.  m  is  odd)  and  the 
“£(2, 2)  extension”  for  even  length  filter.  The  algorithm  consists  of  the  following 
two  steps: 

[1]  extend  the  input  signal  X  to  get  X: 

•  n  is  even  and  m  is  even 

X  =  (^P)  •  •  •  )  ^1)  •  •  •  )  ^ni  •  •  •  j  2:n— p+l)‘ 


•  n  is  odd  and  m  is  even,  if  /  is  a  low-paiss  filter 

X  —  (Xp,  ...,Xi,  Xj,...,  Xn,  X„,  .  .  .  ,  Xn_p) 
and  if  /  is  a  high-pass  filter 

X  —  (Xp,  .  .  .  ,  Xi,  Xj,  .  •  .  ,  X,^,  X,j,  . . .  ,  X„_p^l). 


10 


•  n  is  even  and  m  is  odd,  if  /  is  a  low-pass  filter 

X  —  (Xp4.2,  •  •  •  ,  X2,  Xi,  .  .  . ,  Xti,  X,i_i,  •  •  •  ,  X,i_p)5 

if  /  is  a  high-pass  filter 

X  —  (Xp4.i,  •  •  •  ,  X2,  Xj,  .  .  .  ,  Xn,  X,i_i, .  .  . ,  X71— p— 1)1 

if  /  is  a  low-pass  dual  filter 

X  (Xp-f-l  ,  .  .  .  ,  X2,  Xj,  .  .  .  ,  Xyi,  Xji  — '1^  ,  ,  .  .  ,  X71— p— 1  )  , 

and  if  /  is  a  high-pass  dual  filter 

X  =  (Xp^2?  •  •  •  »  ^2j  j  •  •  •  j  2;^— 1,  •  •  •  j  ^n—p)' 

•  n  is  odd  and  m  is  odd,  if  /  is  a  low-pass  filter 

X  —  (Xp+2»  •  •  •  )  ^2)  ®1)  •  •  •  )  2^71?  ,  •  •  •  ,  X,i_p_i ); 

if  /  is  a  high-pass  filter 

^  —  (Xp+i,  •  •  •  ,  X2,  Xi,  •  •  •  ,  Xre,  Xn— 1,  •  •  •  5  Xn—p); 

if  /  is  a  low-pass  dual  filter 

X  —  (xp^i, . . . ,  X2,  Xi, . . . ,  Xrt,  x,i_i, . . . ,  Xn— p); 

and  if  /  is  a  high-pass  dual  filter 

—  (Xp+25  •  •  •  )  ®2>  ^1»  •  •  •  )  ®n— 1  ?  •  •  •  >  ^n— p— 1  )• 

[2]  apply  the  general  convdown  operator: 

F  =  convdown.general(^,  /). 

4.3  Infinite  (polynomial/zero) 

The  original  signal  Z  is  extended  once  at  the  beginning  of  the  filtering  procedure 
using  a  zero  or  polynomial  of  degree  pdeg  extension.  The  polynomials  are  fit  using 
a  fraction  pf  rac  of  the  data. 

The  infinite  boundary  rule  produces  an  infinite  set  of  wavelet  coefficients.  Only  a 
finite  number  of  coefficients  are  actually  stored.  The  coefficients  which  are  not  stored 
can  be  computed  from  the  stored  coefficients  using  a  zero  or  polynomial  extension. 
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Unlike  the  other  boundary  rules,  the  infinite  boundary  treatment  is  an  “ex¬ 
pansionist”  transform.  This  means  that  in  addition  to  the  desired  Ln/2J  interior 
coefficients  K,  we  also  save  extra  coefficients,  Yt  (left  boundary  coeflBicients)  and  Yr 
(right  boundary  coefficients).  The  number  of  extra  coefficients  depends  only  on  M, 
the  maximum  length  of  analysis  filter  and  synthesis  filter.  For  orthogonal  wavelets, 
the  synthesis  filters  are  the  same  as  the  analysis  filters,  so  M  =  m. 

Note:  In  S-}- WAVELETS,  to  handle  the  boundaxy  coefficients,  the  output  data 
structure  (object  of  cla^s  crystal. list)  is  different  from  other  boundary  rules 
(object  of  cla.ss  crystal. vector). 

Let  p  =  m/2  — 1  and  P  =  M/2  —  1.  Define  d  =  pdegfor  the  polynomial  extension 
and  d  =  — 1  for  the  zero  extension.  The  algorithm  contains  the  following  steps: 

[la]  Initial  polynomial  extension:  if  the  input  signal  is  the  original  data,  proceed 
with  this  step  followed  by  step  2.  Otherwise,  if  the  input  signal  are  the  output 
of  a  previous  convolution  and  down-sampling  operations,  skip  to  step  lb. 

Let  q  =  2(d  P  -f  1)  -|-  p  and  Q  =  max(d  -f  1,  [n  x  pf  rac  +  Ij ). 

•  If  n  is  even,  then 

X  —  (fll,  •  •  •  ,  O5,  •  •  •  )  ^ni  ^1)  •  •  •  5 

•  Otherwise,  if  n  is  odd,  then 

X  —  (fll,  •  •  • ,  flj,  fl-i,  •  •  • )  2:^)  ^1?  •  •  • )  f^9+i) 

The  coefficients  (oi, . . . ,  a,)  are  the  polynomial  or  zero  extrapolations  of  (xi, . . . ,  xq) 
and  the  coefficients  (61,..., 6,)  are  the  polynomial  or  zero  extrapolations  of 
(xn-Q+i, . . . ,  x„).  See  appendix  A  for  details  on  the  polynomial  and  zero  ex¬ 
trapolations. 

Skip  to  step  2. 

[lb]  Extend  X  using  the  existing  boundary  coefficients  (from  previous  filtering 
procedure)  Xe  and  Xr- 

X  =  {Xe,  X,  Xr) 

Let  N  be  the  length  of  X  and  let  q  =  d  +  P  +  1  +  p.  Extend  X  to  obtain  X 
to  obtain 

•  If  A  is  even,  then 

X  —  (fll ,  ...  ,  Qq ,  Xi,  .  ,  .  ,  X  ,  bq ,  .  .  •  ,  bq) . 
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•  Otherwise,  if  N  is  odd,  then 

X  —  (ctl,  •  ■  •  ,  Og,  )  •  •  •  )  ^iVj  ^1)  •  •  •  j  ^g+l)" 

The  coefficients  (ci, . . . ,  a,)  axe  the  polynomial  or  zero  extrapolations  of  (xi, . . . ,  xj+i) 
and  the  coefficients  (6i,...,6,)  cire  the  polynomial  or  zero  extrapolations  of 
. . .  ,Xiv)  (see  appendix  A).  Note  that  pfrac  is  not  used  for  this  extrap¬ 
olation. 

[2]  apply  the  general  convdown  operator  to  X 

Y  =  convdown.general(.A,  /). 

[3]  Let  6  =  d  P  -f  1  and  n'  =  [(n  -I-  1)/2J .  Obtain  the  interior  coefficients  F, 
the  left  boundary  coefficients  Yi,  and  the  right  boundary  coefficients  Yr  by 

y  =  {Vb+x,  •  ■  ■  ,yh+n>) 

Yi  (j/i ,  •  •  • ,  j/ft) 

Yr  —  (yn'+6+lj  •  •  •  )  yn'+2b^ 

On  each  end,  d  +  P  —  I  boundary  coefficients  are  saved. 

4.4  Interval 

An  alternative  to  the  infinite  padding  model  is  to  define  a  slightly  modified  wavelet 
transform  adapted  for  finite  length  signals.  [CDV93]  show  how  to  construct  or- 
thornormai  wavelet  bases  for  intervals  with  compact  support.  In  the  interior  of  the 
interval,  the  basis  functions  are  identical  to  usual  wavelet  basis  functions.  However, 
the  basis  functions  at  the  boundaries  are  given  by  specialized  edge  functions  6^^ 
and  V'J, jt 

The  interval  wavelets  axe  implemented  using  special  filters  at  the  ends  of  the 
series  corresponding  to  the  truncated  and  orthongaiized  wavelet  basis  functions. 

They  are  only  implemented  for  the  discrete  wavelet  transform  with  “symmlets”: 
s4,  s6,  s8,  slO,  sl2,  sl4,  s  16.  The  special  filters  are  stored  as  p-f  1  by  3p  + 2 
boundary  matrices  Be  and  Br-  Sample  size  n  must  divisible  by  2*^  where  J  is  the 
maximum  down-sampling  level  in  the  transform. 

The  algorithm  consists  of  the  following  steps; 

[0]  precondition  X  to  preserve  the  “vanishing  moment”  property  of  the  original 
wavelets.  This  step  is  optional  and  applies  only  to  the  original  series.  Let  Pe 
and  Pr  be  the  left  and  right  precondition  matrices  (transpose  of  S-f  WAVELETS 
functions  left,  precondition  and  right  .precondition  respectively),  then 
(Xi,  .  .  .  ,  Xm/2}  ~  Pti.^li  •  •  •  J  ^m/2)  And  (x,^_,,,i/2+l ,  •  •  ■  ,  Xji)  =  Pe(^X n—m / 2+1  ■>  ■  •  •  1  ^n)  • 
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[1]  extend  X  by  zeros: 


X  =  (0,...,0, 

[2]  apply  the  general  convdown  operator: 

Y  =  convdown.genercLl(X,  /) 

[3]  apply  the  special  filters  to  correct  the  first  p  + 1  and  last  p  + 1  values  of  Y  based 
on  the  first  3p  +  2  and  last  3p  +  2  values  of  X.  Let  Be  and  Br  (S+ WAVELETS 
functions  left .  interval  and  right .  interveil  respectively)  be  the  left  and 
right  boundary  matrices  respectively,  then 


yi  ) 

1  I 

^  Xi  y 

I 

II 

1  Vp+i  J 

^3p+2  / 

/  2/n/2-p+l  ^ 

^  X71— 3p— 1 

1 

=  Br 

yn/2  / 

\  Xn 

See  [CDV93]  for  information  on  obtaining  the  boundary  filters  electronically. 

4.5  Zero  Recursive  Extension 

At  each  filtering  step,  the  coefficients  axe  padded  at  the  beginning  and  the  end  of 
the  signal  with  zeros.  This  boundary  rule  applies  to  arbitrary  sample  size  and  filter 
and  is  implemented  as  a  special  case  of  polynomial  with  pdeg  =  —  1.  The  algorithm 
consists  of  the  following  two  steps: 

[1]  extend  the  input  signal: 

•  n  is  odd  and  /  is  a  low-pass  filter 

p  p+i 

>  X  =  a:i,...,Zn, 

•  otherwise 


A"  =  (0, . . . ,  0,  xi, . . . ,  z„,  0, .  .77^. 
[2]  apply  the  general  convdown  operator: 

Y  =  convdown.general(.A,  /) 
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4.6  Polynomial  Recursive  Extension 

At  each  filtering  step,  the  coefficients  are  padded  at  the  beginning  and  the  end  of 
the  signal  using  a  polynomial  extension  of  degree  d  =  pdeg  fit  to  the  first  d  +  l  and 
Icist  d  +  I  values  respectively.  This  boundary  rule  applies  to  arbitrary  sample  size 
and  filter.  The  algorithm  consists  of  the  following  two  steps: 

[1]  extend  the  input  signal: 

•  n  is  odd  and  /  is  a  low-pass  filter 

X  —  (fli , . . . ,  flp,  , . . . ,  Xfij  ^1 , .  • . ,  ^p+i) 

•  otherwise 

X  ~  (®i)  •  •  •  >  ^1}  •  •  • )  ^1)  •  •  • )  ^p)‘ 

where  (ai,...,ap)  are  the  polynomial  extrapolations  of  (xi, . . .  ,Xd+i),  and 
(bi,... ,  bp)  are  the  polynomial  extrapolations  of  {xn-d,  •  •  •  >  a:„). 

[2]  apply  the  general  convdown  operator: 

Y  =  convdown.general(.A,  /) 

5  Convolution  and  Upsampling  Algorithms 

Let  Y  —  [yi, .  ..^yn)  be  the  input  signal  and  let  /  =  (/i, . . . ,  fm)  be  the  filter.  The 
“generic”  convolution  and  up-sampling  operation  is  given  by 


m 


^2k-i  =  /aj-iyt+fc-i 

»=1 

(5) 

m 

^2k  =  X/  hiVi+k-l 
«=1 

(6) 

Like  the  convolution  and  down-sampling  operator,  the  summation  in  (5)  and  (6) 
needs  to  be  modified  at  the  boundaries. 

Define  the  “convup. general”  operator  by  (5)  and  (6)  with  the  range  of  summation 
restricted  as  follows: 

=  1,2,  ...,n  -p 

where  p  —  m/2  —  1.  The  filtered  vector  is  of  length  2(n  —  p),  which  is  in  general 
smaller  than  the  desired  length  2n  or  2n  —  1.  Therefore,  in  order  to  recover  the 
original  vector,  extra  values  should  be  padded  to  Y  before  calling  convup. general 
operator.  The  extra  values  are  padded  to  the  beginning  and  end  in  the  same  manner 
as  for  the  convolution  and  down-sampling  operator. 
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As  with  the  down-sampling  algorithms  with  the  exception  of  the  reflection  bound¬ 
ary  rule,  we  pad  a  zero  to  /  if  /  is  of  odd  length.  The  filter  padWing  rule  is  as  follows: 
if  /  is  a  low-pass  filter,  /  =  (/,0);  if  /  is  a  high-pass  filter,  /  =  (0,/).  For  the  re¬ 
flection  case,  special  filtering  rules  apply  to  odd  length  filters. 

5.1  Periodic 

When  the  original  series  X  is  assumed  to  be  periodic  with  period  2n,  Y  is  also 
periodic  with  period  n.  Let  p  =  m/2  - 1  and  the  algorithm  consists  of  the  following 
three  steps: 

[1]  extend  Y : 

Y  —  (j/n— p+l  j  •  •  •  1  Vm  yXi  •  •  •  •>  Vni  J/l)  •  •  •  5  Vp) 

[2]  apply  the  general  convup  operator: 

X  =  convup. general(F,  /) 

[3]  select  X  from  X: 

X  —  (Xp+i,  .  .  .  ,  5p+2n)' 

5.2  Reflection 

When  the  filter  is  symmetric/antisymmetric,  Y  is  of  the  same  reflection/periodicity 
property  as  the  original  vector. 

Let  N  be  the  length  of  the  original  vector,  p  =  [m/2j  —  1  and  q  =  [(m  -t-  1)/4J , 
the  algorithm  consists  of  the  following  three  steps: 

[1]  extend  Y  as  follows: 

•  If  iV  is  even,  m  is  even,  then  if  /  is  low-pass  filter 

Y  —  {yqt  Vni  Vm  •  •  •  )  y-n—q-^l ) 

else  if  /  is  a  high-pass  filter 

^  Y  —  {  Vq^  •  •  •  1  y\i  y\i  '  •  '  1  yni  ynt  •  •  •  >  J/n— g+l )• 

#  Else  if  N  is  odd  and  m  is  even,  then  if  /  is  a  low-pass  filter 

Y  —  yi^  •  •  •  ^yny  J/n— 1 )  •  •  •  i  yn—q  ) 

else  if  /  is  a  high-pass  filter 

Y  —  (  Uq,  ,  t/i,  ?/j,...,  yji ,  0,  J/n)  •  •  •  ,  2/n— g+1 )  • 
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•  Else  if  N  is  even  and  m  is  odd,  then  if  /  is  a  low-pass  filter  (not  dual) 

^  (y?)  •  *  •  >  J/l)  •  •  ■  )  Vni  J/n— 1)  •  •  •  >  J/n—p+?— 2)1 

else  if  /  is  a  high-pass  filter  (not  dual) 

^  “  (y^+l)  •  •  ■  ?  1/2}  yii  •  •  • )  Vn}  Vn}  •  •  •  )  J/n— p-1-7— 1)) 

else  if  /  is  a  low-pass  dual  filter 

^  —  (j/g+l  5  •  •  •  j  2/2 )  Vlj  ’  •  •  }  l/n )  J/n  1  •  •  •  >  Vn—p+q—l )  j 

else  if  /  is  a  high-pass  dual  filter  (not  dual) 

y  =  (y„  .  .  .  ,t/l,  I/l,  .  .  .  ,  t/n,  2/n-l,  •  •  •  ,  r/„-p+7-2)- 

•  Else  if  N  is  odd  and  m  is  odd,  then  if  /  is  a  low-pass  filter 

^  {Vq}  •  •  •  ?  2/I5  2/I)  •  •  •  5  Vni  l/m  •  •  •  }  J/n— p+7— 1)) 

else  if  /  is  a  high-pass  filter 

y  —  (2/7+1  )•••  j  y2»  Uli  •  •  •  il/n}  J/n— 1 )  •  •  •  )  J/n— p+g— 1 ) ) 

else  if  /  is  a  low-pass  dual  filter 

y  ~  (j/g+l )  •  •  •  )  2/2?  2/I)  •  •  ■  >  2/n)  Vn—l )  •  •  •  )  yn—p+q—l ) i 

else  if  /  is  a  high-pass  dual  filter 

y  —  (2/7)  •  •  •  )  2/l»  2/l>  •  •  •  )  2/n)  2/n>  •  •  •  »  Vn-p+q-l)- 

[2]  apply  the  general  convup  operator 

X  =  convup.general(y,  /) 
where  /  =  (/,  0)  if  m  is  odd,  and  /  otherwise. 

[3]  select  X  from  X: 

•  if  m  is  multiple  of  4 

X  =  (X2,  .  .  .  ,XjV+l)- 

•  otherwise 

X  =  (xi,  .  .  .  ,XAf). 
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5.3  Infinite  (polynomial/ zero) 

Let  N  be  the  length  of  the  original  vector,  p  =  m/2  —  1  and  P  =  M/2  —  1,  where  M 
is  the  maximum  length  of  analysis  filter  and  synthesis  filter.  The  algorithm  consists 
of  the  following  steps: 

[1]  retract  the  boundary  coefficients  Y(  and  Yr, 

Y  =  (Yi,  Y,  Yr) 

and  apply  the  usual  polynomial  extrapolation  scheme  to  F,  let  qi  =  j^Ez2HSIT£i  j 
and92= 

Y  —  (fli,  • .  • )  ,  y\i  • '  •  t  vni  ^ )  •  •  • )  ^92) 

where  (oi,...,a,i)  are  the  polynomial  extrapolations  of  (yi, . . . ,  and 

(61,..., 6, 2)  are  the  polynomial  extrapolations  of  (j/n-dj •  •  •  jJ/n)-  Note  that 
when  d  =  —  1,  zeros  are  padded  and  when  p  =  P,  91  =  92  =  0,  i.e.  no  padding 
is  needed. 

[2]  apply  the  general  convup  operator, 

X  =  convup.general(F,  /) 

[3]  get  the  interior,  the  left  boundary  and  the  right  boundary  coefficients:  let 
b  =  2{d  +  P  +  1)  —  min(P,p), 

X  —  , . . . , 

X(  =  {xb-d-p,  •  •  •  tib) 

Xr  —  (:^iV+6+l?  ■  •  •  )  ^iV+6+d+J°+l ) 

As  in  convdown  case,  we  save  d  +  P  —  1  boundary  coefficients  on  each  end. 

5.4  Interval 

Let  9  =  3p  +  2.  The  algorithm  consists  of  the  following  steps: 

[\]  apply  the  general  convup  operator  to  Y, 

X  =  convup.genera^F,  /). 

[2]  extend  X  to  the  desired  length  (2n): 


[3]  apply  the  special  filters  to  correct  the  first  q  and  the  last  q  values  of  X  based 
on  the  first  q  and  the  last  q  values  of  Y.  Let  Bi  and  Br  (computed  from 
dwt. matrix)  be  the  left  and  right  boundary  matrices  respectively,  then 


(  y\'\ 

* 

=  Bi 

1 

y<i  J 

^2n-qJr\ 

(  J/n-g+l  ^ 

• 

=  Br 

j 

X2n  yi 

Vn 

1 

The  boundary  matrices  can  be  computed  in  the  following  way:  Let  B  be  the 
transpose  of  the  2q  by  2q  dwt  .matrix  and 


(  B,{L) 

Bc{H)  \ 

V  BAL) 

Br{H)  j 

then  Bt  =  Bt{L)  and  Br  =  Br{L)  if  /  if  a  low-pass  filter;  and  Bt  =  Bi{H) 
and  Br  =  Br{H)  if  /  if  a  high-pass  filter. 


[4]  inverse  precondition  (for  final  step  of  the  synthesis  only). 


5.5  Polynomial/ Zero 


Let  iV  be  the  length  of  the  original  series  X  and  q  =  2max(M/2,  1)  where  M 

is  the  maximum  length  of  analysis  filter  and  synthesis  filter.  The  algorithm  consists 
of  the  following  4  steps: 


[1]  apply  the  general  convup  operator  to  Y : 

X  =  convup.general(y’,  /) 
and  let  n  be  the  length  of  X. 

[2]  extend  X  to  the  desired  length  N, 

•  is  even 

p  p 


•  N  IS  odd  and  /  is  a  low-pass  filter 


X  (0,  .  .  .  ,  0,  X\^  •  *  •  ?  ^n— 2j 
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and  if  /  is  a  high-pass  filter 


p  p+i 

X  =  5 ,  x„,  'o, .  r. ,  0) 

[3]  compute  the  inverse  dwt .  matrix  B  for  the  dual  filter.  The  boundary  matrices 
axe  the  four  comer  matrices  of  B  and  are  different  for  odd  N  and  even  N 
(see  below).  Therefore  there  axe  totally  8  bounadry  matrices,  and  they  axe 
pre-computed  and  stored  in  the  dictionairy. 

[4]  apply  the  left  boundary  matrix  Bt  and  the  right  boundary  matrix  Br  to  correct 
the  boundary  coefficients. 

•  When  N  is  even,  B  is  a  2^  by  2q  matrix.  If  /  is  a  low-pass  filter,  the 
boundary  matrices  Bt  and  Br  are  qhy  q  matrices,  and  defined  by 

Bi  -  1:,] 

Bt  —  -S[(g+1):29, 1:?]; 
axid  if  /  is  a  high-pass  filter, 

Bi  =  5[1:5,  (,+l):2g] 

Bt  =  ■B[(5+i):2g,  (g+l):29]- 


And  then 


/  Xl  \ 


\  J 

^  ^N—q+l  ^ 

\  Xn  J 


=  Bi 


=  Br 


fyi  \ 


^  VN-q+l  ^ 
\  VN  ) 


When  N  is  odd,  5  is  a  2g  -f  1  by  2g  +  1  matrix.  If  /  is  a  low-pass  filter, 
the  boundary  matrices  Bi  and  Br  are  ^-1- 1  by  9  + 1  matrices,  and  defined 
by 

Bi  =  .S[1:(,+1),1:{,+1)] 

Br  =  B((,+i);(2g+l),l;(,+l)]; 

and  then 


^  Xl  ) 

(  y-i-  \ 

* 

=  Bi 

i 

\  ^q-^l  j 

\  ^9+1  / 
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if  /  is  a  high-pass  filter,  Bt  and  Br  are  q  +  l  hy  q  matrices,  and  defined 
by 

Bt  =  -6[1:(,+1),  (,+2):(2,+l)] 

Br  =  -6[(,+1):(2<7+1),  (,+2):(29+1)], 

and  then 


6  Algorithms  for  the  Non-Decimated  Wavelet 
Transform 

Let  X  =  {xi, ,  Xn)  be  the  input  signal  and  let  /  =  (/i, . . . ,  fm)  be  the  filter.  The 
“generic”  convolution  operation  is  given  by 

m 

yk  =  J2  (7) 

«=1 

Since  x,-  is  only  defined  for  i  =  l,...,n,  the  index  k  of  the  output  signal  Y  is 
restricted  to  ^  =  1,2, . . .  ,n  —  m  -t- 1.  Therefore,  the  output  signal  is  always  shorter 
than  the  input  sugnal  X.  In  order  to  obtain  enough  output  coefficients,  m  —  1 
extra  values  should  be  added  to  X.  Currently,  we  have  only  implemented  periodic 
boundary  extension. 


64  Non-Decimated  Wavelets:  Vanilla  Algorithms 

Define  the  convolution  and  dilation  operator  “convdil.general”  by  dilating  filter  / 
and  the  calling  (7).  For  given  input  signal  A  =  (xj, . . . ,  Xn),  filter  /  =  (/i,  •  •  • ,  fm) 
and  a  desired  level  j  (therefore  the  dilation  factor  is  2^),  define  the  dilated  filter 
by  inserting  2^  —  1  zeros  in  each  of  the  adjacent  /,’s,  i.e. 

2-’-l  2J-1 

/2,  •  •  .  fm) 
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and  then  apply  the  “generic”  convolution  operator  (7)  to  X  and 

Let  rrij  be  the  length  of  then  mj  =  (m  —  1)(2-'  —  1)  +  m.  Hence,  for  given 
sample  size  n  aad  filter  length  m,  the  dilation  level  j  is  restricted  to  j  <  log2(^;^). 


6.2  Periodic  Convolution  and  Dilation  Operator 


Forward:  The  algorithm  consists  of  the  following  steps: 
[0]  set  level  j  =  0. 


[1]  at  level  j,  let  p  = 


and  extend  the  input  signal  X 


X  —  (aJn— p+l >  •  •  •  >  ®n)  )  •  •  •  )  ^n>  j  •  •  •  j 


[2]  apply  the  general  convolution  operator  (7): 

Y  =  convdil.generaJ(X,  rev(/),  j). 

[3]  set  j  =  j  +  1.  If  J  <  J,  go  to  step  [1]. 

Inverse:  The  algorithm  consists  of  the  following  steps: 

[0]  set  level  j  =  J. 


[1]  at  level  j,  let  p  =  ,  q  =  and  extend  the  input  signal  X 

X  —  (Xn— p-i-i )  •  •  •  )  ^ni  !  •  •  •  }  Xi,  .  .  .  ,  Ig). 


[2]  apply  the  general  convolution  operator  (7): 

Y  =  convdil.general(X,  /,  j  —  l)/2. 

[3]  set  1.  If  ;  >  0,  go  to  step  [1]. 


6.3  The  A  Trous  Wavelet  Transform 

An  alternative  non-decimated  wavelet  transform  is  given  by  the  “d  Irons"  (“hole”) 
algorithm:  see  [Dut87,  She92,  SMB94].  Like  the  non-decimated  wavelet  transform 
corjiputed  using  nd.dwt,  the  a  irons  algorithm  produces  n  wavelet  coefficients  at 
each  multiresolution  level  for  a  signal  with  n  sample  values.  The  main  difference  is 
that  the  detail  coefficients  in  the  a  irons  algorithm  are  computed  through  simple 
differences  between  the  smooth  coefficients  at  different  levels: 

^j,k  “  ^j—l,k  ~  ^j,k  (^) 

The  detail  coefficients  produced  by  the  nd .  dwt  function  are  computed  using  a  dilated 
high-pass  wavelet  filter. 
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7  Computing  the  DCT 


This  section  gives  discusses  the  algorithms  used  to  compute  the  DCT-II  and  DCT- 
IV  transforms.  For  details  concerning  algorithms  for  the  DCT,  refer  to  the  book  by 
Rao  and  Yip  [RY90]. 

For  a  discrete  signal  /i,  /2,  . . /n,  the  DCT-II  is  defined  as 


2  A,  f{2i-l)kTr' 

cos 


2n 


/:  =  0, 1, . . .  ,n  —  1. 


(9) 


The  scaling  factor  Sk  is  defined  by 


{1  if  k  ^0  and  n 
^  if  =  0  or  n 


(10) 


The  inverse  DCT-II  is  given  by 

/Tn-l 


cos 


n 


k=0 


(2i  —  l)kT^ 
2n 


i  =  1, . . .  ,n. 


(11) 


The  DCT-IV  is  defined  as 


(2i  —  1)(2^  -I-  l)x' 


4n 


^  =  0, 1, . . .  ,n  —  1. 


(12) 


The  inverse  DCT-IV  is  given  by 


n— 1 


IV  ( (2^  ~  1)(2^  -I-  l)7r 


fc=0 


4n 


i  =  1, . . . ,  n. 


(13) 


For  a  discrete  signal  f  =  (/i,  /a?  •  •  •  >  /n),  tbe  DCT’s  and  their  inverses  can  be 
computed  from  the  FFT  of  the  extended  signal  obtained  by  padding  p  zeros  at  the 
end.  Define  the  extension  operator  £^p(x)  for  the  vector  x  =  (xi,  X2»  •  •  ■  >  ^n)  by 


F^p(x) 


Xt  t  =  1,2,  •  •  •  ,n 
0  t  =  n -t- l,n -1- 2,  •  •  •  ,n -f- p 


Define  the  FFT  G{y)  of  a  vector  y  =(yi,  y2,  Vm)  by 

m 

G{y)k  =  Pt  exp  {—i2‘Kk{t  —  l)/m)  A:  =  0, 1,  •  •  • ,  m  —  1  (14) 


£=1 
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7.1  Algorithm  for  the  DCT-II 

For  the  DCT-II,  pad  p  —  n  zeros  and  compute  the  FFT  of  the  signal  axid  take  the 
first  n  frequencies  of  /  =  G{En{f)).  The  DCT-II  coefficients  c  =  (cq,  C2,  . Cn-i), 
are  given  by 

^  f  MO 

1  y!fo  i  =  0  • 

The  inverse  DCT-II  can  be  obtained  by  padding  p  =  3n  zeros  to  the  DCT-II  coef¬ 
ficients  c,  dividing  the  first  DCT  coefficient  cq  by  y/2,  and  taking  the  FFT.  From 
the  FFT  Cfc,  the  original  signal  is  reconstructed  eis  follows: 


fk 


'-Re(c2(k-i)) 

n 


k  =  1, 2,  ...,n. 


(16) 


For  signals  of  length  2'^,  the  DCT-II  is  obtained  from  the  discrete  Haxtley  trans¬ 
form  (DHT).  This  is  a  fast  algorithm  which  avoids  the  need  to  double  the  signal 
length  (as  is  the  case  when  the  DCT-II  is  computed  from  the  FFT)  [Mal86].  The 
discrete  Hartley  transform  [SJBH85]  is  defined  by 


sin 


/  2k{t  —  l)7r' 
\  " 


ft,  (17) 


7,2  Algorithm  for  the  DCT- IV 

For  the  DCT-IV,  pad  p  =  3n  zeros  and  compute  the  FFT  of  the  signal  /  =  G{Esn{^))- 
The  DCT-IV  coefficients  c  =  (cq,  C2,  . . . ,  Cn-i),  are  given  by 


Cfc  =  \/-  cos 
n  \ 


'(2^  +  l)7r\ 
in  J 


Re(f2k+i)  +  sin 


'(2k  -f-  l)7r' 
4n 


Im(f2k+i) 


(18) 


The  inverse  DCT-IV  is  also  given  by  (18),  with  the  coefficients  and  signal  values 
reversed. 


7.3  Boundary  Extension  Rules 

\ 

Since  the  smooth  tapers  extend  beyond  the  ends  of  the  analysis  intervals,  it  is  nec¬ 
essary  to  extend  the  data  at  boundaries.  This  is  similar  to  the  boundary  correction 
algorithms  required  in  wavelet  analysis. 

For  a  taper  of  length  2m  <  n,  it  is  necessary  to  extend  a  signal  of  length  n 
(/ii  /2»  •  •  • )  fn)  by  m  values  on  each  side  to  a  obtain  a  signal  of  length  n  -f  2m 
(/_m+i,  f-m+2,  •••»  fn+m)-  There  are  three  boundary  extension  rules  available  in 
S-fWAVELETS: 
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cp. reflect:  The  signal  is  reflected  at  the  ends  using  the  same  (+,  — )  polarity  as  the 
folding  operator.  The  extended  signal  is 


fi 


f-i+i  for  i  =  -m  +  1,  -m  +  2,  •  •  • ,  0 

<  fi  for  i  =  1, 2,  •  •  • ,  n 

,  /sn-.+l  for  i  =  n  +  1,  n  4-  2,  •  •  • ,  n  +  m 


zero:  The  signal  is  zero  padded,  so  the  extended  signal  is 


fi 


0  for  i  =  —m  +  1,  — m  +  2,  •  •  • ,  0 
<  fi  for  z  =  1,2, •  •  •  ,n 

0  for  i  =  n  +  l,n  +  2, +  m 


periodic:  The  signal  is  eissumed  to  be  periodic  and  wrapped  around  at  the  ends, 
so  the  extended  signal  is 

(  fn-i  for  i  =  -m  +  1,  -m  +  2,  •  •  • ,  0 
fi  =  fi  for  z  =  1,2,  ■••,n 

[  for  z  =  n  +  1,  n  +  2,  ’  •  • ,  n  +  m 

This  is  the  default  boundary  extension  rule  in  S+Wavelets  for  cosine  packet 
analysis. 

Only  the  periodic  boundary  extension  preserves  orthogonality  for  the  boundary 
blocks.  The  boundary  blocks  for  the  cp. reflect  and  zero  extensions  are  only 
nearly  orthogonal. 

Note:  Instead  of  extending  the  data,  Meyer  [Mey93]  suggests  modifying  the 
tapers  at  the  boundary.  By  using  a  discontinuous  “boxcar”  taper,  there  is  no  need 
to  extend  the  signal  beyond  the  boundary.  This  approach  leads  to  an  orthogonal 
transform. 


8  Computing  the  Cosine  Packet  Transform 

Let  X  =  (xi,  ...,Xti)  be  the  signal  and  partition  the  signal  into  p  contiguous  blocks: 
Bi,. . . ,  Bp.  The  jth  block  has  rij  coefficients  and  is  given  by 

=  {fij  5  •  •  ■  J  fij+rij-l) 

where  z'l  =  1  and  ij+i  =  ij  +  nj  iov  j  =  I, . . .  ,p  —  1.  Let  2m  be  the  length  of  the 
taper;  Each  block  Bi  is  of  length  2'^“^  .  Let  2m  be  the  length  of  the  taper  where  m 
is  less  than  or  equal  to  half  of  the  length  of  the  shortest  block.  Let  length  m  vectors 

hint 
hext 


=  \/l  -  iL 
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be  the  weights  of  the  interior  and  exterior  taper  window  respectively.  See  section  B 
below  for  the  tapering  functions  P{t). 

The  cosine  packet  transform  is  computed  by  applying  four  basic  operators  to 
each  block  Bj: 

1.  Use  boundary  rule  to  create  boundary  blocks  Bo  and  Bp+i: 

•  periodic 


Bo  —  Bp 
Bp+i  =  Bo 

Note  that  we  only  need  the  first  m  values. 
•  zero 


m 


•  cp. reflect 


Bo  =  rev(5i) 
Bp+i  =  -rev(5p) 


Note  that  blocks  Bq  and  Bp+i  require  only  m  values. 

2.  Extend  block  Bj  on  the  left  and  right  using  the  neighboring  blocks  Bj_i  and 
Bj+i  to  obtain  the  extended  signal  f  =  (/-m+i,  f-m+2,  fnj+m)- 


/.= 


ki 


for  i  =  —m  -f  1,  —m  +  2,  •  •  • ,  0 

for  i  =  1, 2,  •  •  •  ,nj 

for  i  =  Uj  +  l,nj  +  2,  -  ■  •  ,nj  +  m 


(19) 


For  a  block  of  length  nj,  the  signal  f  is  of  length  2m  +  nj. 


3.  Apply  the  interior  and  exterior  taper  windows  bim  and  6ext  to  block  f  to  obtain 
the  tapered  block  g  as  follows: 

« 


bexti—i  +  !)/»■ 

bmt{i)fi 

fi 

b'mti^^j  ^  "I"  ^')fi 

bext{i  -  nj)fi 


for  i  =  — m  4- 1,  —m  +  2,  ■  •  • ,  0 

for  i  =  1,2,  •  •  •  ,m 

for  i  =  m  +  1,  m  +  2,  •  •  • ,  ny  —  m 

for  i  =  Uj  —  m  +  1,  nj  —  m  +  2,  •••,  Uj 

for  i  =  nj  +  1,  nj  +  2,  •  •  • ,  Uj  +  m 


(20) 
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4.  Fold  m  values  on  each  end  as  follows: 


{ffi-hg-i+i  for  i  =  1,2,- • -,772 

for  i  =  m  +  1,  m  +  2,  •  •  • ,  nj  —  m  (21) 

Qi  ^  ~  771  4"  Ij  Tlj  771  -f-  2,  •  •  •  ,  Tlj 

This  yields  a  “folded”  signal  g  =  (^i,  ^2)  •  •  •  j  Pnj)  of  length  Uj.  This  folding 
is  said  to  have  (+,  — )  polarity,  since  the  reflected  signal  is  added  on  the  left 
boundary  and  subtracted  on  the  right  boundary. 

5.  Apply  the  DCT  to  the  tapered  and  folded  signal  g  to  obtain  the  CPT  coeffi¬ 
cients  Cj  =  Cj,25  •  •  •  j 


9  Computing  the  Inverse  Cosine  Packet  Trans¬ 
form 


The  inverse  cosine  packet  transform  is  obtained  by  reversing  the  steps  of  the  forward 
cosine  packet  transform.  Initialize  the  output  signal  f  to  a  vector  of  n  zeros  (0,  0, 
. . . ,  0).  For  each  block  j  =  1, . . .  ,p,  of  CPT  coefficients  Cj,  do  the  following  four 
steps: 

1.  Apply  the  inverse  DCT  to  the  coefficients  Cj  to  obtain  the  tapered  and  folded 
signal  g. 


2.  Unfold  the  signal  using  the  — )  polarity  to  obtain  the  unfolded  signal  g  of 

length  n  4-  2m. 


9-i+l 


9i  =  i 


9i 


~9nj+{nj-i+l) 


for  i  =  —m  +  1,  — m  +  2,  •  ■  • ,  0 

for  7  =  l,m  -t-  2,  •  •  •  ,nj 

for  i  =  Tlj  +  l,nj  +  2,  ■  ■  ■  ,nj  m 


3.  Untaper  g  to  obtain  the  extended  signal  f: 


7;  = 


+  1) 

9il 

9i 

9i/bmt{nj  -  7  -fl) 

9i  I  77  j) 


for  7  =  —m  +  1 ,  — m  4-  2,  •  •  • ,  0 

for  7  =  1, 2,  •  •  • ,  m 

for  7  =  m  4-  1,  m  4-  2,  •  •  • ,  77j  —  m 

for  7  =  rij  —  m  4- 1, 77j  —  m  4-  2,  •  •  • ,  Uj 

for  7  =  nj  4-  1, 77j  -F  2,  •  •  • ,  77j  4-  777 


(22) 


(23) 


4. 


Separate  f  into  three  blocks:  left  block  Bj^\  central  block  Bj^  and  right  block 


(0,  •  •  •  ,  0,  ,  •  •  •  ,  /o) 
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BJ''^  =  {A>+1)  ■  ■  ■  >0) 

Now  accumlate  the  blocks  as  follows  to  produce  the  blocks  for  the  output  signal: 

Bi^^  B^'^  4"^ 

4'^  4'^ 

_± _ Bj^^  Bj’-) 

Bi  B2  Bp— I  Bp 

The  signal  blocks  Bi,...,Bp  are  constructed  by  adding  the  appropriate  left,  central 
and  right  blocks,  e.g.  B2  =  B^^  +  B2^  +  B^^\ 

Correct  the  boundary  blocks  Bi  and  Bp  based  on  the  boundary  extension  rule: 

•  periodic: 

B^  = 

Bp  =  Bp  +  B[^\ 

•  zero: 

{Bp{n), Bp{n  -  m  +  1)) 

•  cp. reflect: 

{Bril),...,  Brim)) 

.  iBpin),...,Bpin-m  +  l)) 

A  Polynomial  and  Zero  Extrapolation 

In  particular,  for  d  =  0,  ads  are  repeated  values  of  Aveixr,  ■  ■  .,xq).  Note  that  for 
d  —  ~1,  (fll,  ■  •  •  ,  a,q)  (^1)  •  •  •  )  bq)  (O5  •  •  •  )  0)' 


^int(^int  ^ext) 

jBpjn), ...,  Bpjn  -  m  +  1)) 

^int(^int  +  ^ext) 


b> 


int 


iBpin),...,Bpin  -  m  +  1)) 
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For  general  polynomial  extrapolation,  define  a  Vandermonde  matrix  of  (ti, . . . ,  i„): 


/I  1  ...  1  \ 

ti  t2  ...  tn 

V  4.771  ±m  j.m  I 

^2  •  • •  / 


(m-{“l)xn 

and  V  =  Vd{0,  . . . ,  1),  Ve  =  Vd(^, . . . ,  and  K.  =  . 

Then 

/  Oi  \  f  Xi  \ 

=  Vl{VV')-^V 


and 


{h\ 


=  v;{vv')-^v 


\XQ  ) 

(  Xn-Q+\  '' 

V  Xn  / 


<3-1  ) 


B  Tapers  for  Orthogonal  Cosine  Packets 


Very  special  tapering  functions  0j  axe  needed  for  the  cosine  packets  transforms. 
Because  the  cosine  packet  functions  overlap,  only  certain  types  of  tapers  preserve 
orthogonality  in  cosine  packet  analysis  (see  [AWW92,  Wic94]).  S+VVavelets  offers 
7  different  tapers  which  preserve  orthogonality: 


boxcar,  polyl,  poly2,  poly3,  poly4,  polyS,  trig 

The  boxcar  taper  is  discontinuous,  and  is  used  for  the  block  DCT.  The  polyl  is  a 
polynomial  taper  with  one  continuous  derivative.  Likewise,  the  poly2,  . . . ,  polyS 
are  polynomial  tapers  with  2-5  continuous  derivatives.  The  trig  taper  is  based 
on  trigonometric  functions,  and  is  the  smoothest  taper  available  in  S+Wavelets. 
The  default  taper  in  S4- WAVELETS  is  the  poly2  taper. 

A  tapering  function  defined  for  an  interval  I  =  [a,  0]  has  the  form 


m  = 


0 

1 

f^i 

0 


26a  ^ 


P+Sg  —  t  \ 
26^  / 


t  <  a  —  6a 
a  —  6a<t<a  +  6a 
Oi  6a  t  ^  0  —  6p 

0  —  6fi<t<.0-\-6p 

t  0  6(3 


(24) 


where  ^(•)  is  the  left-bell  function  satisfying 


/i(i)  = 


0  t<0 
1  t  >  1 
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and  fi{t)  +  fi{l  —  <)  =  1.  The  parameter  Sa  =  €oA/,  with  Ai  =  ^  —  a,  defines  the 
tapering  region  for  the  left  hand  side  of  the  interval;  and  Sp  =  C/jA/  defines  the 
tapering  region  for  the  right  hand  side,  with  the  constraint  <  1- 

The  specific  form  of  these  tapers  is  follows: 


boxcar  The  boxcar  tapering  function  is  given  by 


0  at  <1/2 
1  if  f  >  1/2 


polyl  ,  . . . , 

by 


polyS  A  polynomial  tapering  function  of  order  p  =  1, . . . ,  5  is  given 


0^ _  i<0 

0<t<l 

1  t  >  1 


where  {bo,  bi,.. . 


,6py  =  ,0y  with 

A.p  —  (uij)o<i,j<p  =  ( 


P  +  j  -1^ 


The  polynomial  tapers  are  available  for  p  =  1,2, ...  ,5,  with  higher  order 
tapers  providing  a  greater  degree  of  smoothness. 


trig  The  tapering  function  for  the  trigonometric  window  is  given  by 


m  = 


0 

sin  (1  —  cos(7rt))) 
1 


t  <  0 
0  <  t  <  1 
t  >  1 


The  trigonometric  taper  is  the  smoothest  taper  available  in  S+ WAVELETS. 
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